Автор: Федоров А.Н.
import sys
import pymol
_stdouterr = sys.stdout, sys.stderr
pymol.finish_launching(['pymol', '-q'])
sys.stdout, sys.stderr = _stdouterr
from pymol import cmd, movie
import time
from IPython.display import Image, HTML
def render(width=640, height=480, filename="/tmp/tmp.png", sleep=1):
cmd.ray(width, height)
cmd.png(filename)
time.sleep(sleep)
return filename
Загрузим все предоставленные в инструкции файлы в текущую папку:
files = ['dppc.itp', 'lipid.itp', 'dppc.gro', 'b.top', 'em.mdp', 'pr.mdp', 'md.mdp']
for f in files:
!wget -q http://kodomo.cmm.msu.ru/~golovin/bilayer/{f}
На основе одного липида созадим ячейку с 64 липидами:
!genconf -f dppc.gro -o b_64.gro -nbox 4 4 4
Посмотрим на наш липид и финальную ячейку в pymol:
%%bash
editconf -f dppc.gro -o dppc.pdb
editconf -f b_64.gro -o b_64.pdb
cmd.reinitialize()
cmd.load('dppc.pdb')
Image(render())
cmd.reinitialize()
cmd.load('b_64.pdb')
Image(render())
Все хорошо, результаты ожидаемы.
В текстовом редакторе в файле b.top установим правильное количество липидов в системе:
// 32 -> 64
Сделаем небольшой отступ в ячейке от липидов, что бы добавить примерно 2500 молекул воды:
!editconf -f b_64.gro -o b_ec -d 0.5
Проведём оптимизацию геометрии системы, что бы удалить "плохие" контакты молекул:
%%bash
grompp -f em -c b_ec -p b -o b_em -maxwarn 2
mdrun -deffnm b_em -v
Начальное значение максимальной силы: Fmax=437970.
Конечное значение максимальной силы: Fmax=616.887.
It works, максимальная сила уменьшилась на порядки.
Добавим в ячейку молекулы воды типа spc:
!genbox -cp b_em -p b -cs spc216 -o b_s
Проведём "утряску" воды:
%%bash
grompp -f pr -c b_s -p b -o b_pr -maxwarn 1
mdrun -deffnm b_pr -v
%%bash
grompp -f em -c b_s -p b -o b_empr -maxwarn 1
mdrun -deffnm b_empr -v
%%bash
grompp -f pr -c b_empr -p b -o b_pr -maxwarn 1
mdrun -deffnm b_pr -v
Посмотрим на b_pr и b_s в pymol:
%%bash
editconf -f b_pr.gro -o b_pr.pdb
editconf -f b_s.gro -o b_s.pdb
cmd.reinitialize()
cmd.load('b_s.pdb')
Image(render())
cmd.reinitialize()
cmd.load('b_pr.pdb')
Image(render())
Нельзя сказать, что здесь заметны какие-то сильные отличия. Все +- похоже, лишь появилось чуть больше молекул воды между липидами.
Для начала протестируем, что все работает локально:
!grompp -f md -c b_pr -p b -o b_md -maxwarn 1
!mdrun -testverlet -deffnm b_md -v -maxh 0.002
Вроде бы все хорошо. Делаем заявку через google-doc.
Path: /home/anfedorov/hw7
tpr file: b_md.tpr
Иии... Ждем.
Подготовим анимацию системы:
!trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20 -pbc mol
!mkdir -p ./frames/
!rm ./frames/*
cmd.reset()
cmd.reinitialize()
cmd.load("b_pbc_1.pdb")
cmd.mpng('./frames/movie')
!ffmpeg -framerate 24 -y -i './frames/movie%04d.png' -c:v libx264 -pix_fmt yuv420p -vf "crop=trunc(iw/2)*2:trunc(ih/2)*2" movie.mp4 > /dev/null
!rm ./frames/*
import base64
with open("movie.mp4", 'rb') as file:
video = base64.b64encode(file.read()).decode()
video = '<video controls alt="test" src="data:video/mp4;base64,{}">'.format(str(video))
HTML(data=video)
Исходная сетка разрушается сразу, бислой начинает формироваться с ~10 модели (t = 4500) и заканчивает к ~25 модели (t = 12000).
!g_traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg
Нормалью к бислою является ось Х. Подсчитаем изменение числа липидов на единицу площади бислоя:
import matplotlib.pyplot as plt
import pandas as pd
df = pd.read_csv('box_1.xvg', header = None, skiprows = 24, sep='\t')
T, X, Y, Z = df[1], df[2], df[3], df[4]
nlipids = 64
plt.plot(T, Y * Z / nlipids)
plt.gca().set(xlabel="Time", ylabel="$\\frac{S_{YZ}}{N_{lipids}}$")
plt.show()
Ожидаемо, число липидов на единицу площади ячейки падает - происходит компактизация системы, липиды организуются в искомый бислой.
!g_sas -f b_md.xtc -s b_md.tpr -o sas_b.xvg -dt 100
Построим зависимость доступной растворителю гидофильной/гидрофобной поверхностей от времени:
df = pd.read_csv('sas_b.xvg', header = None, skiprows = 22, sep='\s+')
T, phobic, philic = df[0], df[1], df[2]
plt.plot(T, phobic, label='Hydrophobic')
plt.plot(T, philic, label='Hydrophilic')
plt.gca().set(xlabel="Time", ylabel="Solvent Accessible Surface (nm\\S2\\N)")
plt.legend()
plt.show()
Теоретически, липидный бислой является более энергетически выгодной структурой, чем просто случайное расположение молекул. Это связано с минимизацией энергетически невыгодных контактов молекул воды и гидрофобных хвостов липидов.
Наши теоретические идеи подтверждаются как качественно(визуализация), так и количественно(график). Видно, что при сборке бислоя резко падает площадь доступной растворителю гидрофобной поверхности. Гидрофобные хвосты липидов группируются вместе, минимизируя число контактов с водой.
Скачаем индекс и запустим gromacs:
!wget http://kodomo.cmm.msu.ru/~golovin/bilayer/sn1.ndx
!g_order -s b_md -f b_md.xtc -o ord_end.xvg -n sn1.ndx -b 45000 -d X
!g_order -s b_md -f b_md.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d X
Строим зависимость:
start = pd.read_csv('ord_start.xvg', header = None, skiprows = 12, sep='\s+')
end = pd.read_csv('ord_end.xvg', header = None, skiprows = 12, sep='\s+')
plt.plot(start[0], start[3], marker='o', label='Start')
plt.plot(end[0], end[3],marker='o', label='End')
plt.gca().set(xlabel="atom", ylabel="oder measure")
plt.legend()
plt.show()
В финальном бислое атомы значительно более упорядочены чем в начале траектории. Это ожидаемо, т.к. бислой - это упорядоченная структура с характерным расположением липидов.
Получилась довольно интересная работа, наглядно показывающая сборку бислоя, фундамента клеточных стенок. Интересно, можно ли смоделировать не только небольшой участок бислоя, но, например, комплекс бислой + мембранный белок и сравнить с готовой pdb моделью (белок, для простоты моделирования, можно брать в уже нативной конформации).
Теоретически, это все еще будет энергетически выгодно, т.к. мембранных белки имеют гидрофобные участки, которые и интегрируются частично/полностью в липидный бислой. А результаты моделирования наглядно бы показали, насколько мы можем доверять МД в подобных задачах.
!jupyter nbconvert hw7.ipynb --to html